%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SVM Decoding
%
% GE_Step1d_clusterstats_ds.m
%
% Statistical evaluation using permutation tests
%
% Modified from GE_Step1d_clusterstats.m
%  by DS 2021-05
%  --Added dashed lines around clusters with p<0.10
%

close all;
clear all;

svm_dir = '/autofs/space/clive_001/users/adriana/GE_SVM';
scripts_dir = sprintf('%s/scripts', svm_dir);
addpath(genpath(scripts_dir));

addpath(genpath('/autofs/space/clive_001/users/adriana/UAG_SVM/scripts/fusionlab_toolbox'));

% Specify the list of subjects:

SubjectNames =  {'GE_05', 'GE_06', 'GE_07', 'GE_08', 'GE_09', 'GE_10', 'GE_11', 'GE_12', 'GE_13', 'GE_14', 'GE_15', 'GE_16', 'GE_17', 'GE_18', 'GE_19', 'GE_20', 'GE_21', 'GE_22', 'GE_23', 'GE_24'}; %'GE_12','GE_13', 
n_subjects = length(SubjectNames);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify the name of the analysis type:

analysis_tag = ''

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify the Regions of Interest ("ROIs", also called "labels":

%roiset_array = {'newSMG', 'newpMTG', 'TempPole', 'FrontalPole', 'CentGyri'};
roiset_array = {'L_CC_1-lh', 'L_cMFG_1-lh', 'L_ParaHip_1-lh', 'L_ParsTri_1-lh', 'L_postCG_1-lh', 'R_ParsOrb_1-rh', 'R_postCG_2-rh', 'R_preCG_1-rh', 'R_preCG_2-rh', 'R_SFG_1-rh', 'L_ITG_1-lh', 'L_ITG_2-lh', 'L_MTG_1-lh', 'L_MTG_2-lh', 'L_ParsOrb_1-lh', 'L_SMG_1-lh', 'L_SPC_1-lh', 'L_STG_2-lh', 'L_STG_1-lh', 'L_TPol_1-lh', 'R_AG_1-rh', 'R_ITG_1-rh', 'R_ITG_2-rh', 'R_ITG_3-rh', 'R_LOC_1-rh', 'R_MTG_1-rh', 'R_MTG_2-rh', 'R_MTG_3-rh', 'R_MTG_4-rh', 'R_postCG_3-rh', 'R_postCG_4-rh', 'R_SFG_2-rh', 'R_STG_1-rh', 'R_STG_3-rh', 'L_SFG_1-lh', 'R_postCG_1-rh', 'R_postCG_5-rh', 'R_STG_2-rh', 'R_SFG_3-rh'};


n_roisets = length(roiset_array);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify the condition pairs:
condition_tag = 'pos'; % 'nvs'|'lvn'|'pos'

switch condition_tag
	case 'nvs'
		%Train all neighbors; test seeds
		categories = {'neighbors'};
		train_conditionA_array = {'11','11','11','11','11','22','22','22','22','33','33','33','44','44','55'};
		train_conditionB_array = {'22','33','44','55','66','33','44','55','66','44','55','66','55','66','66'};
		test_conditionA_array  = { '1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5'};
		test_conditionB_array  = { '2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6'};
	case 'lvn'
		%Separately train lexical and nonword neighbors; test seeds
		categories = {'nonwords'; 'words'}
		train_conditionA_array = {'100','100','100','100','100','200','200','200','200','300','300','300','400','400','500';...
'10','10','10','10','10','20','20','20','20','30','30','30','40','40','50'};
		train_conditionB_array = {'200','300','400','500','600','300','400','500','600','400','500','600','500','600','600';...
'20','30','40','50','60','30','40','50','60','40','50','60','50','60','60'};
		test_conditionA_array = {'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5'};
		test_conditionB_array = {'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6'};
	case 'pos'
		%Separately train based on differential position; test seeds
		categories = {'initial'; 'final'};
		train_conditionA_array = {'1001', '1001', '1001', '1001', '1001', '2001', '2001', '2001', '2001', '3001', '3001', '3001', '4001', '4001', '5001';...
									'1003', '1003', '1003', '1003', '1003', '2003', '2003', '2003', '2003', '3003', '3003', '3003', '4003', '4003', '5003'};
		train_conditionB_array = {'2001', '3001', '4001', '5001', '6001', '3001', '4001', '5001', '6001', '4001', '5001', '6001', '5001', '6001', '6001';...
									'2003', '3003', '4003', '5003', '6003', '3003', '4003', '5003', '6003', '4003', '5003', '6003', '5003', '6003', '6003'};
		test_conditionA_array = {'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5'};
		test_conditionB_array = {'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6'};
    case 'half'
        %Divide data in half randomly (while maintaining proportion of words
        %to nonwords) for training
		categories = {'firsthalf'; 'secondhalf'};
        train_conditionA_array = {'101', '101', '101', '101', '101', '201', '201', '201', '201', '301', '301', '301', '401', '401', '501';...
                                    '102', '102', '102', '102', '102', '202', '202', '202', '202', '302', '302', '302', '402', '402', '502'};
        train_conditionB_array = {'201', '301', '401', '501', '601', '301', '401', '501', '601', '401', '501', '601', '501', '601', '601';...
                                    '202', '302', '402', '502', '602', '302', '402', '502', '602', '402', '502', '602', '502', '602', '602'};
        test_conditionA_array = {'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5'};
		test_conditionB_array = {'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6'};
    case 'third'
        %Divide data in thirds randomly (while maintaining proportion of
        %tokens changed at each position) for training
		categories = {'firstthird'; 'secondthird'; 'thirdthird'};
        train_conditionA_array = {'1010', '1010', '1010', '1010', '1010', '2010', '2010', '2010', '2010', '3010', '3010', '3010', '4010', '4010', '5010';...
									'1020', '1020', '1020', '1020', '1020', '2020', '2020', '2020', '2020', '3020', '3020', '3020', '4020', '4020', '5020';...
									'1030', '1030', '1030', '1030', '1030', '2030', '2030', '2030', '2030', '3030', '3030', '3030', '4030', '4030', '5030'};
		train_conditionB_array = {'2010', '3010', '4010', '5010', '6010', '3010', '4010', '5010', '6010', '4010', '5010', '6010', '5010', '6010', '6010';...
									'2020', '3020', '4020', '5020', '6020', '3020', '4020', '5020', '6020', '4020', '5020', '6020', '5020', '6020', '6020';...
									'2030', '3030', '4030', '5030', '6030', '3030', '4030', '5030', '6030', '4030', '5030', '6030', '5030', '6030', '6030'};
		test_conditionA_array = {'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5';...
									'1','1','1','1','1','2','2','2','2','3','3','3','4','4','5'};
		test_conditionB_array = {'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6';...
									'2','3','4','5','6','3','4','5','6','4','5','6','5','6','6'};
	case 'sop'
		%Position, but reverse training and testing (train on hubs, test on neighbors by position
		categories = {'initial_rev', 'vowel_rev', 'final_rev'};
	case 'biphone'
		categories = {'biphone_first', 'biphone_second'} %  'nonword_first', , 'nonword_second'
        
        test_conditionA_array = {'1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5';...
'1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5'}; %
									
		test_conditionB_array = {'2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6';...
'2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6'}; % 
									
        train_conditionA_array = {'1201', '1201', '1201', '1201', '1201', '2201', '2201', '2201', '2201', '3201', '3201', '3201', '4201', '4201', '5201';...
'1202', '1202', '1202', '1202', '1202', '2202', '2202', '2202', '2202', '3202', '3202', '3202', '4202', '4202', '5202'}; %
									
								
		train_conditionB_array = {'2201', '3201', '4201', '5201', '6201', '3201', '4201', '5201', '6201', '4201', '5201', '6201', '5201', '6201', '6201';... 
'2202', '3202', '4202', '5202', '6202', '3202', '4202', '5202', '6202', '4202', '5202', '6202', '5202', '6202', '6202'}; % 
        
% 		test_conditionA_array = {'1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5';...
% 									'1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5';...
% 									'1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5';...
% 									'1', '1', '1', '1', '1', '2', '2', '2', '2', '3', '3', '3', '4', '4', '5'};
% 		test_conditionB_array = {'2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6';...
% 									'2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6';...
% 									'2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6';...
% 									'2', '3', '4', '5', '6', '3', '4', '5', '6', '4', '5', '6', '5', '6', '6'};
% 		train_conditionA_array = {'1201', '1201', '1201', '1201', '1201', '2201', '2201', '2201', '2201', '3201', '3201', '3201', '4201', '4201', '5201';...
% 									'1301', '1301', '1301', '1301', '1301', '2301', '2301', '2301', '2301', '3301', '3301', '3301', '4301', '4301', '5301';...
% 									'1202', '1202', '1202', '1202', '1202', '2202', '2202', '2202', '2202', '3202', '3202', '3202', '4202', '4202', '5202';...
% 									'1302', '1302', '1302', '1302', '1302', '2302', '2302', '2302', '2302', '3302', '3302', '3302', '4302', '4302', '5302'};
% 		train_conditionB_array = {'2201', '3201', '4201', '5201', '6201', '3201', '4201', '5201', '6201', '4201', '5201', '6201', '5201', '6201', '6201';... 
% 									'2301', '3301', '4301', '5301', '6301', '3301', '4301', '5301', '6301', '4301', '5301', '6301', '5301', '6301', '6301';...
% 									'2202', '3202', '4202', '5202', '6202', '3202', '4202', '5202', '6202', '4202', '5202', '6202', '5202', '6202', '6202';...
% 									'2302', '3302', '4302', '5302', '6302', '3302', '4302', '5302', '6302', '4302', '5302', '6302', '5302', '6302', '6302'};
end

n_condpairs = min(length(train_conditionA_array), length(train_conditionB_array));

% These are for plotting with subplots:
n_hubwords = 6;
n_cols = n_hubwords - 1;
n_rows = n_hubwords - 1;

num_permutation = 1000;

%alpha_array         = [0.05, 0.1];
%cluster_alpha_array = [0.05, 0.1];
alpha_array         = [0.05];
cluster_alpha_array = [0.05];
n_alpha = min(length(alpha_array), length(cluster_alpha_array));

timepoints=1101;

savefigs=1;
%savefigs=0;

zeroed_data_cell=cell(1,n_subjects);
zeroed_ave_data_cell=cell(1,n_subjects);

%accuracy_all = NaN(n_subjects, n_condpairs,timepoints); %store decoding accuracy for all condition pairs for all subjects

time_tag = datetime('now', 'Format', 'yyyyMMdd-HHmm');

for i_category=1:length(categories)

for i_alpha = 1:n_alpha   % try different cluster statistics parameters
    this_alpha = alpha_array(i_alpha);
    this_cluster_alpha = cluster_alpha_array(i_alpha);

  for i_roi=1:n_roisets
    roiset=char(roiset_array(i_roi));
    
    figure('Position',[0, 0, 1650, 1200]);
        
    for i_condpair=1:n_condpairs

        train_conditionA= char(train_conditionA_array(i_category, i_condpair));
        train_conditionB= char(train_conditionB_array(i_category, i_condpair));
        test_conditionA= char(test_conditionA_array(i_category, i_condpair));
        test_conditionB= char(test_conditionB_array(i_category, i_condpair));

        for i_subject=1:n_subjects

            %subject=char(SubjectName(i_subject));
            subject=char(SubjectNames(i_subject));
            subject_dir = char(strcat(subject, analysis_tag));
            input_dir = sprintf('%s/%s/results/', svm_dir, subject_dir);

            results_dir=sprintf('%s/%s%s/results/',svm_dir,subject,analysis_tag);
        
            % Get "accuracy" and "Time" from file:
            inputfile = sprintf('%s/%s_%s_train%sand%s_test%svs%s_Accuracy.mat', input_dir, subject, roiset,train_conditionA,train_conditionB,test_conditionA,test_conditionB);
            load(inputfile);
            if exist('x')
                accuracy = x;
                Time = y;
            end
            zeroed_data_cell{1,i_subject} = accuracy-50;
            
            accuracy_all(i_subject, i_condpair,:) = accuracy; %store decoding accuracy for all condition pairs for all subjects
           
        end %i_subject
        
        % Cluster size analysis:

        %stat = fl_permtestcluster(zeroed_data_cell,'tail','right','statistic','tstat','verbose',0,'numpermutation',1000,'alpha',0.05,'clusteralpha',0.05);
        stat = fl_permtestcluster(zeroed_data_cell,'tail','right','statistic','tstat','verbose',0,'numpermutation',num_permutation,'alpha',this_alpha,'clusteralpha',this_cluster_alpha);
                
        accuracy_ave = mean(accuracy_all(:, i_condpair,:),1);
       
        % plot cluster results
        
        iA = round(str2double(test_conditionA));
        iB = round(str2double(test_conditionB));
        subplot(n_rows, n_cols, (iA-1)*n_cols +  iB-1);

        %plot(Time,accuracy_ave);
        plot(Time,squeeze(accuracy_ave));
        %shadedErrorBar(Time,accuracy_ave_zeroed,accuracy_sem,'patchSaturation',0.16);

        hold on
        %plot(Time,stat.statmap,'LineWidth',5);

        tndx = find(stat.criticalmap);
        %my_color = [.9 .9 .9];
        %my_color = [.3 .3 .3];
        my_color = [.9 .2 .2];
        %plot(Time(tndx),25*ones(size(tndx)),'.','Color',[.7 .7 .7]);
        plot(Time(tndx),75*ones(size(tndx)),'.','Color',my_color);
        %ylim([25,75]);
        ylabel('Accuracy (%)','fontsize',10)
        xlabel('Time (msec)','fontsize',10)
        yline(50);

        titletext=sprintf('Cluster Stats %s: train %s %s, test %s %s',roiset,train_conditionA,train_conditionB,test_conditionA,test_conditionB);
        title(titletext, 'fontsize', 10,'Interpreter','none');
        axis([-200 1100 0 100]);
        legend('off');

    end %i_condpair
    
    %if savefigs
    %    fig_save_stem = sprintf('ClusterStats_%s_condition_pairs_clusteralpha_%s_alpha_%s',roiset,string(cluster_alpha),string(alpha));
    %
    %    results_dir = sprintf('%s/GE_grand_average_results/%s/%s',svm_dir,analysis_tag(2:end),time_tag);
    %    [status, msg, msgID] = mkdir(results_dir);
    %
    %    results_fig = sprintf('%s/%s.fig', results_dir,fig_save_stem);
    %    savefig(results_fig);
    %    results_png = sprintf('%s/%s.png', results_dir,fig_save_stem);
    %    print(gcf,[results_png],'-dpng','-r300');
    %
    %end % if savefigs

    % Cluster Statistics for data averaged over all condition pairs
    % (separately for each subject):
    
    for i_subject=1:n_subjects
        zeroed_ave_data_cell{1,i_subject} = squeeze(mean(accuracy_all(i_subject, :, :),2)) - 50;
    end % i_subject
    accuracy_all_ave = mean(mean(accuracy_all,1),2);
    accuracy_all_std = std(mean(accuracy_all,2),1);
    accuracy_all_sem = accuracy_all_std/sqrt(n_subjects);
    
    stat_ave = fl_permtestcluster(zeroed_ave_data_cell,'tail','right','statistic','tstat','verbose',0,'numpermutation',num_permutation,'alpha',this_alpha,'clusteralpha',this_cluster_alpha);
    
    % Plot the mean accuracy cluster stats in the lower left corner of the figure:
        
    subplot(n_rows, n_cols, (n_rows-1)*n_cols + 1);
    plot(Time,squeeze(accuracy_all_ave));
    shadedErrorBar(Time,accuracy_all_ave,accuracy_all_sem,'patchSaturation',0.16);
    hold on
    %plot(Time,stat_ave.statmap,'LineWidth',5);
    tndx = find(stat_ave.criticalmap);
    %my_color = [.9 .2 .2];
    %plot(Time(tndx),75*ones(size(tndx)),'.','Color',my_color);
    plot(Time(tndx),65*ones(size(tndx)),'.','Color',my_color);
    ylabel('Accuracy (%)','fontsize',10)
    xlabel('Time (msec)','fontsize',10)
    yline(50);

    plot_title = analysis_tag;
    maxlength = 28; % cutoff for too long title information 
    if length(plot_title) > maxlength
        plot_title = plot_title(1:maxlength);
    end
    titletext = sprintf('GE%s: %s (cluster alpha %s, alpha %s)',plot_title, roiset, string(this_cluster_alpha), string(this_alpha));
    %titletext=sprintf('Cluster Stats %s: cluster alpha=%5.3f, alpha=%5.3f',roiset, cluster_alpha, alpha);
    title(titletext, 'fontsize', 10,'Interpreter','none');
    %axis([-200 1100 0 100]);
    axis([-200 1100 25 75]);
    legend('off');

    if savefigs
        fig_save_stem = sprintf('ClusterStats_%s_allcondpairs_clusteralpha_%s_alpha_%s',roiset,string(this_cluster_alpha),string(this_alpha));
        results_dir = sprintf('%s/GE_grand_average_results/%s/%s',svm_dir,categories{i_category},time_tag);
        [status, msg, msgID] = mkdir(results_dir);
        results_fig = sprintf('%s/%s.fig', results_dir,fig_save_stem);
        savefig(results_fig);
        results_png = sprintf('%s/%s.png', results_dir,fig_save_stem);
        print(gcf,[results_png],'-dpng','-r300');
    end % if savefigs
    
    
    % Plot the mean accuracy cluster stats, nicely formatted with gray highlighting:
        
    figure(i_roi*100+i_alpha);
    plot(Time,squeeze(accuracy_all_ave), 'Color', 'b');
    %shadedErrorBar(Time,accuracy_all_ave,accuracy_all_sem,'patchSaturation',0.16);
    %shadedErrorBar(Time,accuracy_all_ave,accuracy_all_std,'patchSaturation',0.16);
    shadedErrorBar(Time,accuracy_all_ave,accuracy_all_std,'lineprops','-b', 'patchSaturation', 0.16);
    
    hold on
    %tndx = find(stat_ave.criticalmap);
    %plot(Time(tndx),65*ones(size(tndx)),'.','Color',my_shading);
    ylabel('Decoding Accuracy (%)','fontsize',10)
    xlabel('Time (msec)','fontsize',10)
    %axis([-200 1100 0 100]);
    %axis([-200 1100 25 75]);
    x_min = -200;
    x_max = 1100;
    y_min = 30;
    y_max = 70;
    axis([x_min x_max y_min y_max]);
    xticks(-100:100:1000);
    yticks(y_min:5:y_max);
    yline(50);
    xline(0);
    legend('off');
    titletext = roiset;
    title(titletext, 'fontsize', 10,'Interpreter','none');

    my_shading = [.5 .5 .5];
	%Get all the clusters, regardless of p-value
%    sig_points = find(stat_ave.statmap>stat_ave.clusteralphamap);
%    clusternum = 1;
%    all_clusters{1} = sig_points(1);
%    for j_sig = 2:length(sig_points)
%        if sig_points(j_sig) == sig_points(j_sig - 1) + 1
%            all_clusters{clusternum} = [all_clusters{clusternum}; sig_points(j_sig) ];
%        else
%			clustersize(clusternum) = length(all_clusters{clusternum});            
%			clusternum = clusternum + 1;
%            all_clusters{clusternum} = sig_points(j_sig);
%        end
%    end

%	clustersize(clusternum) = length(all_clusters{clusternum});
%	pvalues = sum(stat_ave.clustermassdist>clustersize)./1000;
%
%    for i_clust = 1:length(all_clusters)
 %       iic = all_clusters{i_clust};
%		if pvalues(i_clust) < 0.1        
%			if cell2mat(cellfun(@(x) any(cellfun(@(y) isequal(x,y), stat_ave.clusters)), all_clusters(i_clust), 'uni', false)) %if the cluster is in the significant clusters cell array
%        	    harea = area([Time(iic(1)) Time(iic(end))], [100, 100], 'FaceColor', my_shading, 'LineStyle', 'none');
%        	else
%        	    harea = area([Time(iic(1)) Time(iic(end))], [100 100], 'FaceColor', 'none', 'LineStyle', '--', 'EdgeColor', my_shading);
 %       	end
%		else
%			continue
%		end
%        alpha(0.2);
%    end %i_clust
%    hold off;

% Replaced by above code in order to also input p<0.10 cluster highlighting
    for i_clust = 1:length(stat_ave.clusters)
        iic = cell2mat(stat_ave.clusters(i_clust));
        harea = area([Time(iic(1)) Time(iic(end))], [100 100], 'FaceColor', my_shading, 'LineStyle', 'none');
        alpha(0.2);
    end % i_clust
   hold off;
    
    if savefigs
        fig_save_stem = sprintf('ClusterStats_%s_clusteralpha_%s_alpha_%s',roiset,string(this_cluster_alpha),string(this_alpha));
        results_dir = sprintf('%s/GE_grand_average_results/%s/%s',svm_dir,categories{i_category},time_tag);
        [status, msg, msgID] = mkdir(results_dir);
        results_fig = sprintf('%s/%s.fig', results_dir,fig_save_stem);
        savefig(results_fig);
        results_png = sprintf('%s/%s.png', results_dir,fig_save_stem);
        print(gcf,[results_png],'-dpng','-r300');
        stat_filename = sprintf('%s/cluster_stats_%s', results_dir, roiset_array{i_roi});
        save(stat_filename, 'stat_ave');
    end % if savefigs

  end %i_roi

end %i_alpha

end %i_category

